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Abstract 

We present a rigorous procedure to derive coarse-grained red blood cell (RBC) models, which lead 
to accurate mechanical properties of realistic RBCs. Based on a semi-analytic theory linear and non- 
linear elastic properties of the RBC membrane can be matched with those obtained in optical tweezers 
stretching experiments. In addition, we develop a nearly stress-free model which avoids a number of pitfalls 
of existing RBC models, such as non-biconcave equilibrium shape and dependence of RBC mechanical 
properties on the triangulation quality. The proposed RBC model is suitable for use in many existing 
numerical methods, such as Lattice Boltzmann, Multiparticle Collision Dynamics, Immersed Boundary, 
etc. 

Key words: atomistic modeling, dissipative particle dynamics, spectrin model 



1 Introduction 

Recent experiments on the red blood cell (RBC) to probe its mechanical properties include micropipette 
aspiration ([D, 0) and RBC deformation by optical tweezers 13, 0, 0)- These experiments provide clear 
evidence that RBCs subject to deformations are characterized by a non-linear mechanical response. The 
healthy human RBC assumes a biconcave shape with an average diameter of 7.8 \im. RBCs have relatively 
simple structure (@, 13) comprising of a membrane filled with a liquid cytosol of fixed volume. The RBC 
membrane consists of a lipid bilayer with an attached cytoskeleton formed by a spectrin protein network 
linked by short actin filaments. The lipid bilayer can be considered nearly viscous and area-incompressible 
<H), 

while the attached spectrin network is mainly responsible for the membrane elastic response providing 
RBC integrity as it undergoes severe deformations in narrow capillaries as small as 3 \im in diameter. This 
relatively simple RBC structure can be considered as an excellent system to study and rigorously model its 
complex behavior and response. To this end, a number of numerical models have been developed recently, 



including continuum descriptions da |<j HfJ, Hi ) and discrete approximations at the spectrin molecular level 



(12, 13) as well as at the mesoscopic scale (14, 15, 0, Fully continuum (fluid and solid) type of 
modeling often suffers from difficulties in coupling nonlinear solid motions and fluid flow and from ex- 



cessive computational expense. Therefore, "semi-continuum" modeling dlfj, lilt) of deformable particles 
is developing rapidly and is typically based on the immersed boundary or front-tracking techniques. Here 
a membrane is represented by a set of points which move in a Lagrangian fashion and are coupled to an 
Eulerian discretization of fluid domain. Continuum models often do not consider membrane thermal fluc- 
tuations present at the mesoscopic and microscopic scales. On the other hand, detailed spectrin molecular 
modeling of RBCs is limited due to extreme computational demands. Therefore, in this work we will focus 
on accurate mesoscopic modeling of RBCs. 

There exist several mesoscopic methods \vh for modeling deformable particles and in 

particular RBCs. Dzwinel et al. (16) model RBC as a volume of elastic material which has an inner 
skeleton. This model does not take into account the main structure concept of RBC (a membrane filled with 
a fluid), and therefore it cannot capture the proper dynamics, for example the tumbling and tank-treading 
behavior in shear flow (fl8lfl9|) . The other three aforementioned methods (fi^[l5ll"l7^ employ a very similar 
approach where the RBC is represented by a network of springs in combination with bending rigidity and 
constraints for surface-area and volume conservation. Dupin et al. (15J) couple the discrete RBC to a 
fluid described by the Lattice Boltzmann method |2fj). Their results are very promising, however the main 
disadvantage of their model is the fact that thermal fluctuations are not considered, which are of importance 
in RBC rheology and dynamics. Noguchi and Gompper (14) employed Multiparticle Collision Dynamics 



(1211) and presented encouraging results on vesicles and RBCs. Pivkin and Karniadakis (1171) used Dissipative 
Particle Dynamics (DPD) (|22|) for a coarse-grained RBC model which we will use as the starting point 
of our work. Specifically, we will develop a generalized RBC model with major improvements in its 
mechanical properties. 

The paper is organized as follows. In the next section we present details of the RBC model. Section [3] 
provides a semi-analytical theory of the RBC membrane elastic properties and compares the RBC stretch- 



General coarse-grained red blood cell models: I. Mechanics 



2 



ing deformation with experimental data. We conclude in section 0] with a brief discussion and detailed 
suggestions for model development. 



2 Red blood cell model 



The membrane model structure is analogous to the works in (1141 115L 1171) . It is defined as a set of points 
with Cartesian coordinates {x^}, i G l...N v , that are vertices in a two-dimensional triangulated network on 
the RBC surface. The vertices are connected by N s edges represented by springs, which form Nt triangles. 
The free energy of the system is given by 

^({xj}) — Vin— plane ^bending ~i~ V area + V vo l ume . (1) 

The in-plane free energy term includes the energy of springs, U s , and may contain the elastic energy 
stored in the membrane as follows 

Vin— plane = ^ ] U s (lj) + ^ ] ^4^"' 
jei...N s feel— JV t k 

where lj is the length of the spring j, is the area of the k-th triangle, and the constant C q and exponent q 
should be properly selected. Different spring models can be used here and we will discuss performance of 
some of them in section[3] However, we would like to highlight two nonlinear spring models: the wormlike 
chain (WLC) and the finitely extensible nonlinear elastic (FENE) spring, the attractive potentials of which 
are given, respectively, by 

k j^T^Ij jymix 2^27 kg o r o-i 

UWLC = "j : , UpENE = -TT l max lo g L 1 ~ X ' ( 3 ) 

4p 1 — x I 

where x = l/l max £ (0, 1), l m ax is the maximum spring extension, p is the persistence length and k s is 
the FENE spring constant. Note that when the distance between two connected points approaches l max , 
the corresponding spring force goes to infinity, and therefore limits the maximum extension to l m ax- It is 
important to point out that both WLC and FENE springs exert purely attractive forces, thus they produce 
a triangle area compression, while the second term in equation (O provides a triangle area expansion. The 
stressless (or minimum energy) state corresponds to an equilibrium spring length Iq and depends on the 
introduced parameters. The relation among these parameters and the equilibrium length can be derived by 
an energy minimization argument (1121) or by setting the Cauchy stress obtained from the virial theorem to 
zero (|23T) . We obtained the following expressions for WLC and FENE springs, respectively, 



-,wlc = V3A q +1 k B T(4x 2 -9x + 6) FENE = V3A q +1 k s 

/q 4pql max (l-x )2 ' * " q(l-xl) 



where xq = Iq/ l m ax and Aq = \/3Zq /4. These formulas allow us to calculate the strength of the second 
term in equation Q for the given equilibrium length and spring parameters. Another choice is to consider 
a spring with a specific equilibrium length (e.g., harmonic spring, WLC or FENE in combination with 
a repulsive potential), and then set C q to zero. We now introduce a repulsive force defined as a power 
function (POW) of the separation distance I as follows 

fpow(l) = ^, m > 0, (5) 

where k p is the force coefficient and m is the exponent. The combination of WLC or FENE with POW 
defines a spring with an equillibrium length, and will be called WLC-POW and FENE-POW, respectively. 
The strength k p can be expressed in terms of the equillibrium length Iq and the WLC or FENE parameters 
by equating the corresponding forces. The combination of WLC or FENE with the in-plane energy in 
equation (0 will be denoted as WLC-C and FENE-C throughout the paper. 
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The bending energy is defined as 

Vbending = V h [1 - C0s(9j - 9 )] , (6) 



J61...AT S 



where fe& is the bending constant, 9j is the instantaneous angle between two adjacent triangles having the 
common edge j, and 9q is the spontaneous angle. 
The area and volume conservation constraints are 

_ k*(A - A^f v kgjAj - A f 

Varea- tt + ^ 2A ' 1 ' 

jei...N t u 

K{v-_yt? 

A 

where k a , k^ and k v are the global area, local area and volume constraint constants, respectively. The 
terms A and V are the total area and volume of RBC, while Aff* and Vq° 1 are the desired total area and 
volume, respectively. Note, that the above expressions define the global area and volume constraints while 
the second term in equation (7a) corresponds to local area dilatation. 

In order to obtain the forces corresponding to the above energies we use the following formula 

fi = -dV({xi})/dxi, i £ l...N v . (8) 

Exact force expressions can be derived analytically from the defined energies, however for brevity we do 
not present them in this paper. 



3 Mechanical properties 

Mechanical properties of RBCs were measured in a number of experiments by micropipette aspiration 
(CIS) and RBC deformation by optical tweezers 

aaa. 

The reported shear modulus /j,q lies between 2 and 
15 fiN/m and the bending modulus k is between 1 x 10~ 19 and 7x 10~ 19 J, which corresponds to the range 
of 23 — 163 ksT based on the normal body temperature T = 36.6°C. In addition to some uncertainties in 
the experiments, the discrepancies in the measurements arise in part from applying simplified geometrical 
models to extract values from the measured forces as the precise geometry is often not known. In such 
cases, accurate numerical modeling can provide a valuable aid in experimental parameter quantification. 

In recent optical tweezers stretching experiments (0,11) the RBC behavior was modeled using a hyper- 
elastic material model and the finite element method (FEM). From the FEM simulations the range for the 
membrane shear modulus of hq = 5 — 12 fiN/m was obtained. This corresponds to the Young's modulus 
of Y = = 15 — 36 ixN/m due to the three-dimensional membrane model. Dao et al. (|23r) performed 
coarse-grained molecular dynamics (CGMD) simulations of the spectrin-level cytoskeleton which yielded 
a worse approximation to the experimental stretching response in comparison with FEM. They derived 
the first-order approximation of the shear modulus (j,q and the area-compression modulus K for a regular 
hexagonal network of springs expressed through spring parameters. Eventhough the shear modulus in the 



FEM and CGMD simulations was matched, it is clear from figure 8 of (1231) that FEM and CGMD sys- 
tems have different Young's modulus as the slopes in the linear elastic deformation regime are different. 
In addition, their estimated area-compressibility modulus was K = 2// which yields the Poisson's ratio 
of v = 1/3, while the membrane was nearly incompressible. However, for an incompressible material 
one can find that u = 1 for a two-dimensional membrane model and K — > oo. We have confirmed that 
their analytical results are correct but they appear to be incomplete because not all model contributions are 
considered for the membrane elastic properties estimation, which can explain the inconsistency found. We 
will explain this further in the next section. 
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3.1 Linear elastic properties 

Our starting point is the linear analysis of a two-dimensional sheet of springs built with equilateral triangles 
as described in detail in (I23I) . Figure Q] (left) shows an element of the equilateral triangulation with vertex 
v placed at the origin. The stress for the area element S (from the virial theorem) is given by 



T ^ = ~2A 



a a ap H T—Oabp H 



(b a - a a )(bp - dp) 



' A q+l 



at 



A tot 



N t A) k d (A - A) 



An 



lap, 



(9) 



where /(•) is the spring force, a, (3 can be x or y, N t is the total number of triangles and A\ ot = NAq. In 
general, Nt cancels out and the global and local area contributions to the stress can be combined together as 
— (k a + kd)(Ao — A)/A()5 a p. Note, that the linear analysis in (23) didnot take into account the global and 
local area contributions to the stress which significantly affects the final results. The linear shear modulus 
can be derived by applying a small engineering shear strain 7 to the configuration in figure Q] (left) and 
taking the first derivative of shear stress /in = 7pp| 7 =o- The shear deformation is area-preserving, and 
therefore only spring forces contribute to the membrane shear modulus. For different spring models, we 
obtained the following expressions for n$: 



,,WLC~C 



V3k B T 

^plmaxXQ \4(1 
y/3k s 



3 A xq 

- + 4:X + — ~ 

4 2(1 
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(10a) 

(10b) 
(10c) 
(lOd) 



The linear elastic area compression modulus K can be calculated from the area expansion with the 
resulting pressure given by 



1 



ij~XX ~t~ Tyy) 



3Z/(Q 
AA 



+ q 



C q , (k a + k d )(A - A) 



AH 



+ 



(11) 



The compression modulus K is defined as 

dP 



K 



d log (A) 



1 dP 



A=A 



2 d log (I) 



1 dP 



l=lo 



2 d log (x) 



(12) 



X=XQ 



Using equations (fTTT> and (fl2l we derive the linear area compression modulus for different spring models 
as follows 



K 



WLC-C 



V3k B T 



q + 



K 



1 

2 

1 — x, 
K WLC-POW _ 



(Ax 2 - 9x + 6) + 



FENE-C 



q+ 1 + 



1 



1 + 2(1-X ) 
1 — Xq 



+ k a + kd, 



31 



+ & a + & d 



A" 



FENE-POW 



2fi WLC-POW + ka+kdj 
2fl FENE-POW + K + ^ 



(13a) 

(13b) 

(13c) 
(13d) 
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Note, that if q = 1 we obtain the expressions K WLC ~ C = 2p^ LC ~ c + k a + k d and K FENE ~ C = 
2fi FENE ~ c + k a + kd. Generally, for a nearly incompressible sheet of springs the area constraint coeffi- 
cients have to be large such that k a + ^» 1, and thus K S> /j,q. 

The Young's modulus Y for the two-dimensional sheet can be expressed through the shear and area 
compression modulus as follows 



and the Poisson's ratio v is given by 

K - no 



1, if K ->oo. (15) 



The above expressions are consistent with the incompressibility assumption enforced through the condition 
k a + kd S> 1. In practice, we use the value of k a + k d = 5000 which provides a nearly incompressible 
membrane with the Young's modulus of about 2% smaller than the asymptotic value of 4/io (^o = 100). 
All the analytical expressions for p,Q, K and Y were numerically verified by shearing, area expanding 
and stretching experiments of the regular two-dimensional sheet of springs. In addition, it is important to 
note that the modeled sheet appears to be isotropic for small shear and stretch deformations, however it is 
anisotropic at large deformations. 

3.2 Membrane bending properties 

In this section we discuss the correspondence of our bending model to the macroscopic model of Helfrich 



(|24r) given by 

E = 7T / (Ci + C2 - 2C fdA + k 9 I dC 2 dA, (16) 

2 J A J A 

where C\ and C2 are the local principal curvatures, Co is the spontaneous curvature, and k c and k g are the 
bending rigidities. 

We base the derivation on the spherical shell. Figure[T](right) shows two equilateral triangles with sides 
a, the vertices of which rest on the surface of a sphere of radius R. The angle between their normals ni and 
ri2 is equal to 6. For the spherical shell we can derive from equation ( fT6l ) E = 87rA; c (l — Co / C\) 2 +A-Kk g = 
87rA; c (l — R/Rq) 2 + Airkg, where C\ = C2 = l/R and Co = l/i?o- For the triangulated sphere we 
have E t = N s kb[l — cos(9 — 9q)] in the defined notations. We expand cos{9 — 0q) in Taylor series 
around (6 — 9q) to obtain E t = N s kb(9 — 9q) 2 /2 + 0((9 — #o) 4 )> where we neglect high-order terms. 
From figure Q] (right) we find that 2r w 9R or 9 = and analogously #0 = ^ R . Furthermore, 

A sphe re = ^R 2 fn N t A = = and thus a 2 /R 2 = 8ttV3/N s . Finally, we obtain E t = 

Nshi^ - 7M- ) 2 / 2 = ^m^i 1 ~ R / R o) 2 = f7f U - R/Ro) 2 - Equating the macroscopic bending 
energy E for k g = —Ak c /3, Co = (25J) and E t gives us the relation kj, = 2k c /y/3, which is the same as 
derived in the continuum limit in (1251) . The spontaneous angle #0 is set according to the total number of 
vertices N v on the sphere. It can be shown that cos (6) = 1 — -1/4) = — Wtt) / (V3N S — 6ir), 

while N s = 2N V — 4. The corresponding bending stiffness k^ and the spontaneous angle #0 are given by 

h = T^ 9 ° = C0S [ WN v -2)-to )- (17) 

3.3 RBC triangulation 

The average unstressed shape of a single red blood cell measured in the experiments in (@) is biconcave 
and is described by 



4(x 2 + y 1 

J5 2 



x 2 + y 2 (x 2 + y 2 ) 2 



clq + ai — —I-, 1- 02 



D 2 D 4 



(18) 
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where D = 7.82 fim is the cell diameter, a = 0.05179025, oi = 2.002558, and a 2 = -4.491048. The 
area and volume of this RBC is equal to 135 /j,m 2 and 94 \iw?, respectively. We have investigated three 
types of triangulation strategies: 

• Point charges: N v points are randomly distributed on the sphere surface, and the electrostatic problem 
of point charges is solved while the point movements are constrained on the sphere. After steady state 
is reached the sphere surface is triangulated and conformed to the RBC shape according to equation 

OB. 



Gridgen: the RBC shape is imported into commercially available grid generation software Gridgen 



(26) which performs the advancing front method for the RBC surface triangulation. 

• Energy relaxation: First, the RBC shape is triangulated following the point charges or Gridgen meth- 
ods. Subsequently, the relaxation of the free energy of the RBC model is performed while the vertices 
are restricted to move on the biconcave shape in equation (fl~8l) . The relaxation procedure includes 
only in-plane and bending energy components and is done by flipping between the two diagonals of 
two adjacent triangles. 

The triangulation quality can be characterized by two distributions: (i) distribution of the link (edge) length, 
(ii) distribution of the vertex degrees (number of links in the vertex junction). The former is characterized 
by the value d{l) = cr(l)/T, where I is the average length of all edges, and a (I) is the standard deviation. 
The latter defines the regularity of triangulation by providing the relative percentage of degree-re vertices 
re = l...n max - Note that the regular network, for which the mechanical properties were derived, has only 
degree-6 vertices. Table Q] presents the mesh quality data on average for different triangulation methods. 
The better mesh quality corresponds to a combination of smaller d(l), higher percentage of degree-6, and 
smaller percentage of any other degree vertices, and is achieved for larger number of points N v . At this 
point, it seems that the best quality is reached with the free energy relaxation method while the worst is the 
Gridgen (advancing front method) triangulation, which will be discussed further below. 

3.4 Dissipative Particle Dynamics modeling and scaling to real units 

We will model RBC with Dissipative Particle Dynamics (DPD), a mesoscale method used for simulations 

j I 

of complex fluids and soft matter, see (22) for details. We now outline the scaling procedure which relates 
DPD non-dimensional units to real units. First, we choose the equilibrium spring length Iq = iff in DPD 
units, and the superscript D denotes "DPD" and [iff] = r c , where r c defines DPD length scale. Another 
parameter we are free to select is the imposed shear modulus fj,Q = fiff with [nff] = ^— = ^ kB y , which 
will provide a scaling base. Use of WLC and FENE springs requires to set the maximum extension length 
lm ax , however it is more convenient to set the ratio r mu it = iff /I max- Further in the paper we will show 
that the choice of r mu i t does not affect the linear elastic deformation, but it governs the RBC non-linear 
response at large deformation. For given iff, [iff and r mu i t we can calculate the required spring parameters 
for a chosen model using equations (lOa-d). Then, the area-compression modulus K D and the Young's 
modulus Y D are found for the calculated spring parameters and given area constraint parameters (k a and 
fcrf) using equations (13a-d.[T4"l). At this point, we can define the length scaling based on the cell diameter 
Dff = (Lff + L*>)/2, where [Dff] = r c and L x , L y are the cell diameters in x and y directions found from 
the equilibrium simulation of a single cell using the previously obtained model parameters. The length 
scaling based on iff appears to be inappropriate, because, in general, the cell dimensions would depend on 
the relative volume to area ratio and to some extent on the present triangulation artifacts (discussed later 
in text). As an example, we can define RBC and a spherical vesicle with the same iff, while the cell sizes 
would greatly differ. However, in general, Dff is proportional to iff for fixed volume to area ratio. The real 
RBC has the average diameter Dff = 7.82 fim (superscript R denotes "real"), and therefore the following 
length scaling is adapted 

D R 

r c = -^[m}. (19) 
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Due to the fact that we will perform RBC stretching simulations, it is natural to involve the Young's 
modulus into the scaling as the main parameter. Matching the real and model Young's modulus Y D ( fca -p 

_ yR ( fc ^) provides us with the energy unit scaling as follows 

yR r 2 yR / tjR\ 2 

(kBT) D = - T ,-3,(k B T) R = - f , ( ° ) (k B T) R . (20) 



Y D m 2V ° ' Y D \ Dq 

After we determined the DPD energy unit (as an example for the human body temperature of T = 36.6°C), 
we can calculate the bending rigidity in DPD using the energy unit and equation (fT71) . In addition, we define 
the force scaling by 

r c Y D Dq m v ' 

Note that for the stretching simulations here we do not need to explicitly define mass and time scaling as 
we are not interested in stretching dynamics. 

3.5 RBC stretching: success and problems 

We perform RBC stretching simulations and compare results with the experimental data of RBC defor- 
mation by optical tweezers (5). Here, we assume that the real RBC has diameter D R = 7.82 \im. The 
aforementioned FEM simulations of RBC membrane ( 5) showed an agreement with the experimental data 
for p, R = 5.3 fiN/m, however we find that a slightly better correspondence of the results is achieved for 
/j, R = 6.3 jiN/m and Y R = 18.9 fxN/m, which we select to be the targeted properties. Table [2] shows 
a set of the required RBC parameters. The triangulation for all N v was performed using the free energy 
relaxation method. The imposed Young's modulus for all cases is Y D = 392.453, which is about 2% 
lower than that in the incompressible limit Y D = 4 / Uq ) = 400. Using equation d20l ) we find the energy 
unit (hsT)® based on {kBT) R at the normal body temperature T = 36.6°C. The bending rigidity k c is 
set to 2.4 x 10 -19 J, which seems to be a widely accepted value and is equal to approximately 56(kBT) R . 
The total RBC area A ot is equal to ^^(Zq 5 ) 2 , where Nt is the total number of triangle plaquettes with 
the area A = ^{Iq) 2 - Note that for all triangulations used in this paper N t = 2N V - 4. The total RBC 
volume V^ ot is found according to the following scaling V^ ot /(A^) 3 / 2 = V R /(A R ) 3 / 2 , where V R = 94 
and A R = 135 fim 2 according to the average RBC shape described by equation (fl~8l) . 

The modeled RBC is suspended in a solvent which consists of free DPD particles with number density 
n = 3. Note that macroscopic solvent properties (e.g., viscosity) are not important here, because we 
are interested in the final cell deformation for every constant stretching force. Thus, we allow enough 
time for the RBC to reach its final deformation state without close monitoring of the stretching dynamics. 
Meanwhile, the solvent maintains the temperature at the constant value of {ksT) . 

Figure[2]shows a sketch of the red blood cell before and after deformation. The total stretching force F R 
is in the range 0...200 pN, which can be scaled into DPD units F S D according to equation (l2Tb . The total 
force F<P is applied to N + = eN v vertices (drawn as small black spheres in figure 0) of the membrane with 
the largest x-coordinates in the positive x-direction, and correspondingly — F S D is exerted on N- = N + 
vertices with the smallest x-coordinates in the negative x-direction. Therefore, a vertex in N + or iV_ is 
subject to the force = ±F^ /(eN v ). The vertex fraction e is equal to 0.02 corresponding to a contact 
diameter of the attached silica bead d c = 2 \im used in experiments. The contact diameter was measured 
as ^maxjj \yf — + maxjj \y~ — yj\j /2, where yf, y^ and y^ , yj are the y-coordinates of vertices 
in N + and iV_, respectively. The simulations for the given force range were performed as follows: (i) 
M = 16 is chosen, which defines the force increment AF R = 200 pN/M with corresponding AF/ 3 . (ii) 
The loop % = 1...M is run with the stretching force % • AF^ during time 2r each. The time r is long 
enough in order for RBC to converge to the final stretching state for the given force. Thus, the time [0, r] 
is the transient time for convergence, and during time [r, 2r] the deformation response is calculated. The 
axial diameter Da is computed over time r as \x max — x m i n \, where x max is the maximum x position 
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among the N+ vertices, while x m in is the minimum among N-. The transverse diameter Dt is calculated 
as 2 x max i=1 jv„ \J {Ui — c y ) 2 + (z~i — c z ) 2 , where c y , c z are the y and z center of mass coordinates. 

Figure [3] presents the RBC stretching response for different number of vertices N v (left) and spring 
models (right) with RBC parameters from tabled also included are experimental results (5) and the coarse- 
grained RBC model results of (1171) . We find an excellent agreement of the simulation results with the 
experiment independently of the number of vertices or spring model. A noticeable disagreement in the 
transverse diameter may be partially due to the experimental errors rising from the fact that the optical 
shape measurements were performed from a single observation angle. RBCs subjected to stretching may 
rotate in y-z plane which was noticed in numerical simulations, and therefore measurements done from a 
single observation angle may result in underprediction of the maximum transverse diameter. However, the 
simulation results remain within the experimental error bars. The solid line in figure [3] corresponds to the 
coarse-grained RBC (17) of similar type employing the WLC-C model. In (fl7b the derivation of linear 
elastic properties did not include a contribution of the area constraint, which results in Young's modulus 
underprediction on the order of 50%. From the region of small near-linear deformation (0 — 50 pN) it 
is clear that the solid line corresponds to a membrane with a larger Young's modulus compared to the 
experiment. In addition, in order to compensate for the error in the estimated membrane elastic properties 
the ratio r mu i t was set to 3.17, which results in near-linear elastic deformation, and ignores a non-linear 
RBC response at large deformations. Finally, it is worth commenting that the FENE-C model appears to 
be less stable (requires a smaller time step) at large deformation due to a more rapid spring hardening 
compared to WLC-C. The WLC-POW model performs similar to WLC-C, however a weak local area 
conservation (kd > 0) may be required for stability at large deformations as it mimics the second in-plane 
force term in equation © for the WLC-C model. 

Despite the demonstrated success of the RBC models, several problems are remaining due to a not 
stress-free membrane. Figure [4] shows the RBC response for different triangulation methods: free energy 
relaxation triangulated WLC-C N v = 500 RBC stretched along lines with distinct orientation angles (left) 
and RBC response for models with different triangulation (right). While RBC triangulated through free 
energy relaxation method gives satisfactory results with difference in the stretching response on the order 
of 5 — 8%, RBCs triangulated by other methods show much greater discrepancy with the experiment. In 
addition to that, RBCs triangulated by point charges and Gridgen methods require to set the bending rigid- 
ity to 300(/cbT) r and 200(/cbT) r , respectively, in order to maintain the equilibrium biconcave shape, 
while the bending rigidity of the real RBC is about 56(/cbT)- r . Here, a lower bending rigidity may result 
in relaxation to the stomatocyte (cup) shape (13). Moreover, figure 0] shows that these models have higher 
effective elastic modulus than that predicted as they are subject to a higher membrane stress at equilib- 
rium due to triangulation artifacts. Also, they appear to give a stronger stretching anisotropy (10 — 15%) 
compared to the free energy relaxation method. 



3.6 Stress-free membrane model 

To eliminate the aforementioned membrane stress artifacts we propose a simple modification to the de- 
scribed model. For each spring we define F i = l...N s set to the spring lengths after the shape triangula- 
tion, since we assume it to be the RBC equilibrium state. We define accordingly l l max = Iq x r mu u and A 3 G 
j = l...N t for each triangular plaquette. The total RBC area A ot = Ylj=i N t an ^ tne tota l volume 
V tot is calculated from the RBC triangulation. Then, we define the average spring length Tq = ^i=i n s 
and the average maximum spring extension as l max = Tq x r mu u, which are used in the linear elastic 
properties estimation using equations (14c,d and 17c,d). Here, we omit the WLC-C and FENE-C models 
because it may not be possible to define a single in-plane area expanding potential (the second force term 
in equation (|2])) for a triangle with distinct sides. However, for the WLC-POW and FENE-POW models the 
individual equilibrium spring length can be simply defined. Thus, the WLC and FENE spring parameters 
(ksT/p and k s ) will be the same for all springs calculated through Tq, T max and iXq, while the power force 
coefficient k p in equation ([5]> will be adjusted in order to set the given equilibrium spring length. 

We perform tests using the WLC-POW model for different triangulation methods and number of ver- 
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tices. Figure [5] presents simulation results for N v = 500 with different triangulations (left) and a range 
of the number of vertices N v from 100 to 27344 (right). A substantial improvement is observed when 
compared with the results in figure [4] (right). Note that the stress-free model, when probed along different 
stretching directions results in deviation in the stretching response on the order of 1% for the free energy 
triangulation method and about 3 — 5% for the other triangulation techniques. The stretching response 
for different number of vertices gives an excellent agreement of the results with the experiment. Here, 
N v = 27344 corresponds to a spectrin-level of RBC modeling as in (1131) . while N v = 100 — 500 is 
highly coarse-grained RBC. Eventhough the coarse-grained model of N v = 100 yields correct mechanical 
deformation results, it may not provide an accurate smooth RBC shape description which can be of impor- 
tance for the dynamics. We propose that the minimum N v to be used for the RBC model should be about 
250 - 300. 

Dependence of the RBC deformation response on the ratio r mu i t and on the number of vertices N + , 
N- (figure [2]) is shown in figure [6] As we mentioned before small RBC deformations are independent of 
the ratio r mu n, however at large deformation this parameter plays a significant role and governs the non- 
linear RBC response. In addition, figure [6] (right) shows that the RBC response is sensitive to the fraction 
of vertices (shown in percent) to which the stretching force is applied. It is equivalent to changing d c in 
figure [2 which characterizes the attachment area of silica bead in the experiments. 



3.7 Comparison with a single spectrin tetramer 

It is rather remarkable that RBCs can be accurately modeled with just a few hundred points, which is about 
hundred times computationally cheaper than the spectrin-level RBC model, where N v ~ 27000 (23). At 
the spectrin-level of RBC modeling, each spring represents a single spectrin tetramer, and therefore the 
spring force WLC-POW should mimic the spectrin tetramer deformation response. We are not aware of 
any experimental single spectrin stretching results, however this has been done by means of numerical 
simulation using a coarse-grained molecular dynamics (CG-MD) approach in (1271) . Figure [7] presents a 
single spectrin tetramer stress-strain response compared with the spring force of the spectrin-level RBC 
model. We find a remarkable agreement of the results. Here, we assume that the maximum extension 
spring length is 200 nm as in the CG-MD simulations of (|27l) . This corresponds to Iq = 91 nm with 
fmuit = 2.2 from where an appropriate length and force scaling can be calculated as done in the previous 
sections. 



4 Summary 

We presented general coarse-grained RBC models represented by a network of springs in combination with 
bending rigidity, area, and volume conservation constraints. The modeled RBC accurately captures elastic 
response at small and large deformations, and agrees very well with the experiments of the RBC stretch- 
ing with optical tweezers. The linear elastic properties of the RBC membrane are derived analytically, 
and therefore no manual adjustment of the model parameters through numerical tests is required. We also 
proposed a stress-free RBC model which leads to triangulation-independent membrane properties, while 
the conventional RBC model suffers from non-zero local stresses which result in triangulation dependent 
deformation response and equilibrium shape. The model was tested for different levels of coarse-graining 
starting from the spectrin-level modeling (about N v = 27000 vertices) and ending with only N v = 100 
vertices for the full membrane representation. However, we suggest that the minimum number of vertices 
to be used for the RBC membrane should be about N v = 250 — 300, while the lower N v may not accu- 
rately represent RBC smooth shape which is of importance for the RBC dynamics. In addition, we found an 
excellent agreement of single spring force in case of the spectrin-level model with the spectrin tetramer re- 
sponse obtained from the coarse-grained molecular dynamics simulations. The proposed model is general 
enough, and therefore can be easily applied in many numerical methods, such as semi-continuum meth- 
ods (Immersed Boundary and Advanced Front Tracking), mesoscopic methods (Lattice Boltzmann and 
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Brownian Dynamics), and mesoscopic particle methods (Dissipative Particle Dynamics and Multiparticle 
Collision Dynamics). 

Here, we summarize the procedure for the RBC model. First, we obtain triangulation of the equilibrium 
RBC shape defined by equation (|T8T > for the given number of vertices N v . This triangulation sets the 
required equilibrium lengths for the springs, triangle areas and the total RBC area and volume. Second, 
we choose the modeled membrane shear modulus /i , and area and volume constraint coefficients (eq. 
(7a,b)). This defines our RBC model parameters using equations (lOa-d, 13a-d, 14 and 17) with the average 
equilibrium spring length, and will set a scaling to the real units using equations (20,21). In addition, we 
need to define the length scaling (eq. (19)) based on the RBC diameter. We suggest to obtain the RBC 
diameter through an equilibrium simulation rather than assuming it from the analytical RBC shape (eq. 
(Q~8|>) as they may be slightly different depending on the relative contributions of in-plane elasticity and 
membrane bending rigidity. After these two simple steps, the linear elastic properties of the model will 
match those of the real RBC. In addition, we need to mention that in the case of strong RBC deformations 
we may need to adjust the spring maximum extension length which governs the non-linear RBC response. 
However, it is convenient to set the ratio r mu i t = lo/lmax = 2.2 for the WLC springs and r mu i t = 2.05 for 
the FENE springs. We emphasize that the described procedure does not involve any parameter adjustments 
through a number of numerical tests. 

The spectrin stretching comparison provides additional justification of using the spring model for ac- 
curate RBC deformation response. From these results we can draw the conclusion that: an appropriate 
spring model for RBC should have the maximum allowed extension length, in the neighborhood of which 
the spring force rapidly hardens in order to prevent further membrane strain. In view of this, the harmonic 
spring used in (114) gives an adequate response at small deformations but it will not capture a non-linear 
RBC deformations. Furthermore, the neo-Hookean spring used in (|15[) provides a good RBC stretching 
response but it may fail at very large deformations. At this point, an experimental confirmation of the 
single spectrin tetramer stress-strain relation would be of great interest. 

This paper is the first part of our work on accurate coarse-grained RBC modeling. This simple coarse- 
grained model is inexpensive and it accurately captures the RBC mechanical properties. The second part of 
our work concerns RBC rheological properties and dynamics. It will demonstrate the importance of treating 
the RBC membrane as a viscoelastic material, and will present studies of the complex RBC dynamics and 
rheology. 
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5 Tables 



Method 


d(l) 


degree- 6 


degree-5 and degree-7 


other degrees 


point charges 


[0.15,0.18] 


90% - 95% 


5% - 10% 


0% 


Gridgen 


[0.13,0.16] 


45% - 60% 


37% - 47% 


3% - 8% 


free energy relaxation 


[0.05,0.08] 


75% - 90% 


10% - 25% 


0% 



Table 1: Mesh quality for different triangulation methods. 



Model 


N v 


l D 
L o 




D D 
u q 


r c 


Tmult 


k a 


k d 


k v 


#0 


q or m 


WLC-C 


500 


0.56 


100 


8.267 


1.0 


2.2 


5000 





5000 


6.958° 


1 


WLC-C 


1000 


0.4 


100 


8.285 


1.0 


2.2 


5000 





5000 


4.9° 


1 


WLC-C 


3000 


0.23 


100 


8.064 


1.0 


2.2 


5000 





5000 


2.821° 


1 


FENE-C 


500 


0.56 


100 


8.265 


1.0 


2.05 


5000 





5000 


6.958° 


1 


WLC-POW 


500 


0.56 


100 


8.25 


1.0 


2.2 


4900 


100 


5000 


6.958° 


2 



Table 2: RBC parameters. 



6 Figure Legends 



Figure 1: An element of the equilateral triangulation (left) and two equilateral triangles placed on the 
surface of a sphere of radius R (right). 

Figure 2: RBC sketch before and after deformation. 

Figure 3: Computational results for different N v (left) and spring models (right) compared with the 
experiments m <S) and the coarse-grained (CG) RBC model in (1171) . 

Figure 4: RBC stretching along lines with different orientation angles (left) and triangulation methods 
(right) compared with the experiments in qI). 

Figure 5: Stress-free RBC model for different triangulation methods with N v = 500 (left) and number 
of vertices with the energy relaxation triangulation (right) compared with the experiments in (Q). 

Figure 6: The stretching response of the stress-free RBC model for different ratio r mu i t (left) and num- 
ber of vertices in percents which are subject to the stretching force (right) compared with the experiments 
in<@). 



Figure 7: A single spectrin tetramer stress-strain response (27) versus the spring force of the spectrin- 
level RBC model. 



7 Figures 
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